## ---------------------------------------------------------------------- #
## Figure 7
## Christopher J. Fariss
## "Yes, Human Rights Practices Are Improving Over Time"
## American Political Science Review
## ---------------------------------------------------------------------- #
rm(list = ls())

year <- 1946:2015


NAMES <- c(
"+ CIRI: Disappearance (4 items)",
"+ PTS: Amnesty (5 items)",
"+ CIRI: Extrajudicial Killing (6 items)",
"+ Hathaway: Torture (7 items)",
"+ CIRI: Torture (8 items)" ,
"+ PTS: State (9 items)"
)



for(i in 1:6){
    
    MODEL <- i
    
    
    if(MODEL==1) load("progressive_no_events_summary_disap.RData")
    if(MODEL==2) load("progressive_no_events_summary_amnesty.RData")
    if(MODEL==3) load("progressive_no_events_summary_kill.RData")
    if(MODEL==4) load("progressive_no_events_summary_hathaway.RData")
    if(MODEL==5) load("progressive_no_events_summary_tort.RData")
    if(MODEL==6) load("progressive_no_events_summary_state.RData")
    
    
    if(MODEL==1) test <- read.csv("progressive_no_events_summary_amnesty.csv")
    if(MODEL==2) test <- read.csv("progressive_no_events_summary_disap.csv")
    if(MODEL==3) test <- read.csv("progressive_no_events_summary_kill.csv")
    if(MODEL==4) test <- read.csv("progressive_no_events_summary_hathaway.csv")
    if(MODEL==5) test <- read.csv("progressive_no_events_summary_tort.csv")
    if(MODEL==6) test <- read.csv("progressive_no_events_summary_state.csv")
        
    
    tmp <- subset(test, grepl("theta", test$X), select=c(mean, sd))
    
    
    data  <- data[order(data$id),]
    
    data$latentmean <- tmp$mean
    data$latentsd <- tmp$sd
    
    year <- YEAR <- 1946:2015
    
    temp1 <- temp2 <- NA
    for(i in 1:length(year)){
        m <- data$latentmean[data$YEAR==year[i]]
        v <- data$latentsd[data$YEAR==year[i]]
        temp <- apply(matrix(rnorm(matrix(1,ncol=1000, nrow=length(m)), mean=m, sd=v),ncol=1000, nrow=length(m)), 2, mean)
        temp1[i] <- mean(temp)
        temp2[i] <- sd(temp)
    }
    latentmean <- ifelse(year<1976, -999, temp1)
    #latentmean <- temp1
    latentsd <- temp2
    
    data <- as.data.frame(cbind(YEAR, latentmean, latentsd))
    
    COLOR <- c(2,3)
    COLOR <- c(1,1)
    
    temp <- matrix(-999, nrow=5000, ncol=length(year))
    
    # take 5000 draws from a normal distribution for each year of the data, which fills 5000 rows in each i column of the data
    for(i in year-1945){
        temp[,i]  <- rnorm(5000, mean=data$latentmean[data$YEAR==year[i]], sd=data$latentsd[data$YEAR==year[i]])
    }
    
    #par(mar=c(2,3,0.2,2), cex=1.2, font=2, font.lab=2)
    par(mar=c(2,5,2,2), cex=1.2, font=1, font.lab=1)
    
    # calculate the mean value of each column in the matrix; i.e., the mean of each column of 5000 draws
    point.estimates <- apply(temp, 2, mean)
    
    # plot the country-year point.estimates
    plot(point.estimates, ylim=c(-0.6,1.20), yaxt="n", xaxt="n", pch=19, col=grey(.25), xlab="", ylab="", cex=.8, type="n", xlim=c(1,70))
    axis(1, at=(1:70), label=rep("", length(1:70)))
    axis(1, at=(seq(5,65,10)), label=seq(1950,2010,10), lwd=1.5)
    axis(2, at=c(-0.75,-0.50,-.25,0,0.25,0.50,0.75, 1.0), las=2)
    #axis(2, at=c(-1,-0.5,0,0.5,1.0), las=2)
    mtext(side=2, "Latent Physical Integrity", line=3.5, cex=1.7)
    mtext(side=4, "More Abuse", font=1, line=0.3, cex=1.7, at=-0.15)
    mtext(side=4, "More Respect", font=1, line=0.3, cex=1.7, at=0.65)
    arrows(75.25, -.50, 75.25, -.65, xpd = TRUE, lwd=3, length=.15)
    arrows(75.25, 1.05, 75.25, 1.20, xpd = TRUE, lwd=3, length=.15)
    
    abline(h=-.25,lty=2, col=grey(.75), lwd=2)
    abline(h=0,lty=2, col=grey(.75), lwd=2)
    abline(h=.25,lty=2, col=grey(.75), lwd=2)
    abline(h=.5,lty=2, col=grey(.75), lwd=2)
    abline(h=.75,lty=2, col=grey(.75), lwd=2)
    abline(h=1.00,lty=2, col=grey(.75), lwd=2)
    
    # add 95% credible intervales to the plot
    for(i in 1:length(year)){
        lines(c(i,i), c(apply(temp, 2, quantile, .025)[i], apply(temp, 2, quantile, .975)[i]), lwd=2, col=grey(.6))
    }
    points(apply(temp, 2, mean), col=grey(.6), pch=21, bg=grey(.95), cex=.8)
    
    
    
    
    mtext(side=3, NAMES[MODEL], line=0.3, cex=1.7)
}









